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1 . INTRODUCTION 


1.1. The Problem 

The NASA Lewis Laboratory is conducting an investigation 
into the influence of zero and reduced gravitational accel- 
eration on diffusion flames, with a view to improving under- 
standing of fires in space vehicles. Experiments have been 
conducted by NASA Lewis personnel on the burning in air of a 
jet of methane emerging from a • cylindrical tube; the air is at 
rest relative to the tube, apart from motion induced by the 
jet; the whole apparatus can be dropped down a tower, so that 
zero-gravity operation may be investigated. 

In experiments in which the apparatus is first at rest, and 

then suddenly dropped, the flame length passes from one 

(nearly) steady value to another value, after passing through 

smaller transient values. The re-establishment of a steady 

state lasts about 0.5 seconds. The tube radius is in the 

range 0.05 to 0.44 cm, and the flow rate of methane is in the 
3 

range 1 to 8 cm per sec. The maximum velocity at the tube- 
exit axis is around 1200 cm per sec. The whole experiment lasts 
for about 5 seconds. The pressure is atmospheric, and the air- 
and gas- supply temperatures are around room temperature. The 
flow is laminar. 

The NASA Lewis Laboratory wishes to compare its experiment- 
al findings with numerical solutions of the differential equa- 
tions believed to govern the processes, with a view to confirm- 
ing ..qualitative understanding of the process, to determine 
the quantitative accuracy of numerical predictions, and to 
establish a mathematical model of the process'- for subsequent 
use as a predictive and exploratory tool. 

During recent years CHAM Ltd has developed a range of analytical 
models and numerical techniques to study fluid-dynamic , heat 



-6 - 

discussion of the analytical results obtained in this study. 
The final chapter is devoted to concluding remarks concern- 
ing the computational accuracy, and efficiency of the comput- 
er code, together with an assessment of the physical real- 
ism of the results. In each case, suggestions are made con- 
cerning means by which they may be improved. 

1.2. Previous Work 

i — — » 

The earlier analytical study of steady laminar diffusion flames 
can be categorised in terms of firstly, semi-empirical app- 
roaches based upon correlations of specific' data (see typ- 
ically Refs 1 and 2) and secondly, special solutions of the 
conservation equations, (e.g. Refs 3 and 4). These have been 
summarised in References 5 and 6. These approaches either ig- 
nored the effects of gravity. or provided inadequate agree- 
ment with experiment when gravitional effects are significant. 

More recent developments using more precise numerical tech- 
niques allied with greater computing power in modern comput- 
ers have enabled more sophisticated analytical models to be 
studied. These developments recognise the need to solve the 
equations for mass, momentum and energy in circumstances 
which take into account interdependence of the velocity, 
temperature and concentration fields. These methods also use 
finite-difference forms of the governing differential equa- 
tions which reduces the restrictions on the generality of 
the boundary conditions and provide a framework for rapid 
development of the analytical models to include for 
example the effects of gravity. Edelman et al.(Ref 7) dev- 
eloped a model with sufficient generality that further 
insights could be gained concerning the effects of gravity 
on steady laminar diffusion "flames . Good agreement was 
achieved with experimental results in normal-g flames, how- 
ever, predictions were not so good under conditions of Hero 
gravity. 
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transfer and combustion problems. These models and tech- 
niques have been applied to a variety of problems, both 
steady and unsteady, with and without chemical reaction. 

The major purpose of the study described in this report 
was to demonstrate that these modelling and numerical tech- 
niques could be used to study the transient flame problem 
which NASA’s Lewis Laboratory has been studying experiment- 
ally, and, furthermore, to demonstrate that a computer code 
based upon the analytical model and numerical techniques 
can be developed to give good quantitative as well as qualit- 
ative agreement with experiment. 

CHAM used’ its EASI (Elliptic Axi^ymmetric Integrator) code 
as a starting point. The specific derivative of this code 
which was assembled for this project described herein was 
termed EGOMAC (Effect of Gravity on Methane-Air Combustion) 

The program of work has been done in two parts. Eirstly EGOMAC' 
was assembled and run for a particular choice of geometry 
and fuel-supply conditions. When the code was seen to 'produce 
qualitatively correct results for the transient flame, the 
second part of the work was devoted to a systematic investig- 
ation of the effects on the computed solutions of:- computat- 
ional grids; time steps; - outer boundary conditions; number of 
iterations; relaxation factors; and solution order of equat- 
ions, The code was additionally exercised for two other fuel 
.flow rates, and two other tube diameters. In addition the 
sensitivity of the solutions to changes in the specified physi- 
cal chemical constants, was also studied. 

The remaining two parts of the introductory chapter are dev- 
oted to a description of previous work, and to outlining the 
specific new features embodied in the CHAM approach to this 
problem. The second chapter describes the mathematical anal- 
ysis upon which the computer code is based. The third and 
fourth chapters contain respectively a presentation and 
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Personnel at CHAM Ltd have' been developing numerical models 
and solution algorithms for predicting axially flowing 
reacting flows. Patankar and Spalding (Reference 8) describe 
the GENMIX code which has been widely used for diffusion 
flames under ’'parabolic" conditions, i.e. those of steady 
flow without upstream propagation of influences from down- 
stream. 

Edelman et al, concluded that axial diffusion of mass, energy 
and momentum which introduced "elliptic" effects (i.e. 
propagation upstream of influences from downstream) becomes 
--important at low-Reynolds-Number , zero-gravity conditions. 
This indicates an explicit requirement of the analytical 
method solving the transient laminar flame problem, viz that 
it should be an elliptic method. 

CHAM has developed a basic elliptic code (EASI) which is 
-^ell suited to this problem, and has therefore used it as 
the basis for the development of the EGOMAC code which is 
described in this report. 

1.3. Outline of the Present Contribution 

The EGOMAC code contains all of the computational advantages 
contained in the basic EASI code i.e. rapid convergence and 
accuracy without the excessive computing times usually assoc- 
iated with elliptic codes. The main contribution of the work 
presented in this report^ is to demonstrate that the EGOMAC 
code will efficiently compute the hydrodynamic and thermo- 
dynamic fields associated with the laminar flame under trans- 
ient gravity conditions. This demonstration, renders the 
validity of the novel features embodied into the code self- 
evident. These novel features include :- 

1. Inclusion of axial diffusion in the physical model which 
in certain cases is the mechanism which keeps the flame 
alight. 

Provision of lateral pressure gradients in the model. 


2 . 
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3. The use of hybrid differencing scheme to enhance both 
accuracy and stability. 

4. The use of implicit rather thah explicit procedures for 
calculating all variables in the flow field (this is 
the heart of the SIMPLE Algorithm). 

During the study described herein, it should be noted that 
additional computational efficiency was achieved by ensuring 
that the initial conditions, corresponding to a lg steady 
flame were computed using the GENMIX code. In this way, it 
was possible to ensure that the transient calculation started 
from a fully converged and precise definition of all relev- 
ant physical quantities. This together with the significance 
of the novel features contained in the EGOMAC code will be 
described in more detail in later chapters. 
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2. MATHEMATICAL ANALYSIS 

2 . 1 Differential Equations and Boundary Conditions 
Differential Equations 


Equations are solved for the following eight dependent 
variables : 

o the components of velocity in the axial and radial 
(x and r) directions - u and v, 

© the static pressure - p, 

e ' the enthalpy h; - 

© the mixture fraction - f - defined as the proportion 
of fluid present at any point which originated, from 
the nozzle, 

© the mass fraction of fuel (methane) - rn^ u , 

© and the radiation-flux quantities - and R^, which 

are defined in terms of the radiative heat fluxes in 
the positive and negative x direction (K and L) and the 
positive and negative r direction (I and J), as follows : 

R x = l (K + L) (2.1-11) 

R t = $‘ (I + J) (2.1-2) 

The conservation equations for most of the above quantities 
(viz. for all except p, R x and R r , which are considered later) 
are solved in the general form appropriate to an unsteady, 
axisymmetrie flow given below : 

9 (V<jA + P (Hm + L <L Fr( - n, <LtY] = S, 

it r 1 * Sr 4 
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In this- equation , t is the time, $ is the general dependent 
variable and p is the density. 

The diffusion coefficient for <j> (T^) is evaluated as follows, 

Tor velocities, F^ = y (2.1*4) 

and, for other variables, F^ = y/c^ (2.1—5) 

Here, y stands for the laminar viscosity of the fluid 

(which is calculated as described in Section 2.2 below), 

and 0 ^ is the Prandtl or Schmidt number. In the present 

work, a-. , a. and cr are taken to be constant and equal to 
h 1 m fu 

0.6, except that for one calculation (case 26) is set to 
0.7. 


Values are given in Table 2-a for the source terms S^. The 
symbol g there stands for the gravitational acceleration in 
the negative x direction (see Fig. 2-a). In the source term 
forh, a- is the absorption coefficient, and E is the black- 
body emissive power, which is related to T (the local 
absolute fluid temperature) and the Stef an-Boltzmann 
constant (cr), by : 


E = aT 1 * 


( 2 . 1 - 6 ) 
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s + 

u 

-a+£W+£rrf r i‘?l)-»<’ 

TT 

5r r w ' Br y r 1 

' \ 


2a(ft x + R r )--4-aEI 


O 

• m fa 

2. "f m fu m oX exp(-A/&T) 


' Table 2-a : List of Source Terms 
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Finally, in the source term for m fu (which is a consequence 
of presuming that the reaction-rate is given by an Arrhenius 
expression) and A are constants whose values are set as 
described in Section 3. while m is the concentration of 0 9 . 

The pressure p, is deduced, by the method described in 
Section 2.3, from the mass-continuity relation, which : 


It + + ll (f^r) = o 

it fx - r 3r 


(2.1-7) 


Finally, the equations for R x and R y are 


A ( r„ iB-) + 5. = o 

dx - ax ' 


(2.1-8) 


ii r R jju + s r = o 

. I Jr 1 r dr j 


(2.1-9) 


where. 


r. _ & +• s) 


( 2 . 1 - 10 ) 


P R “ *■/( « + s + i/r) 


( 2 . 1 - 11 ) 


S. = a E 4- 5 CR X 4 + <2.1-12) 
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S Rip = aE -f s (R x + R r )/l - (q + s)R 1> (2.1-13) 

The symbol s here stands for the scattering coefficient. 

Eqns. (2.1-8) and (2.1-9), it should be observed, can be 
regarded as especially simple forms of the general <f> 
equation (Eqn. (2.1-3); they lack the time derivative (since 
radiation is not stored), the convection terms (since 
radiation is not conveyed by bulk movement), and the diffusion 
term in the direction normal to that of the radiation flux. 

Boundary Conditions 


The solution domain is shown in Fig. 2-a. Boundary conditions 
are applied at all four boundaries as described below. 

At the nozzle, the rate of mass inflow is specified; the 
axial-velocity u profile is parabolic, and the radial 
velocity (v) is zero. The mixture fraction f is unity by 
definition; and when, as in the case here, the incoming 
fluid is pure methane, iri Ju is also unity. The enthalpy li is 
calculated from the specified fluid temperature as described 
in Section 2.21. 

At the remainder of the lower (AB) boundary, and at the outer 
(BC) boundary, h is deduced from the specified ambient 
temperature; u, f and mj- u are zero at both boundaries, and 
the lateral velocity v is zero at the lower boundary. At 
the outer (BC) boundary, the v velocities are set so as to 
account for the lateral mass inflow by entrainment. 

The distribution of mass inflow was previously determined 
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from solutions with the boundary- layer procedure, GENMIX 
(Ref. 8), and is tabulated below. However, it is important 
to remark that, for the calculations reported, the outer 
boundary was chosen to be sufficiently remote from the 
flame, that the effect of the conditions prescribed there 
on solutions within the flame was negligible.* 


x cm. 

_ 1 

kg . s 

0.00 

0 

0.40 

1.64^10“ 6 

0.90 

4 . 10 x 10~ 6 

2.10 

9.56 x 10" 6 

3.18 

1;44*10~ 5 

4.46 

1 . 95* 10~ 5 

5.08 

2. 18 x 10~ 5 

6.49 

2.65* 10~ 5 

7.27 

2.92 x 10~ 5 

8.26 

3.24*10“ s 

, 9.40 

3.57*10~' s 

20. '70 

3.96 x 10~ s 

12.20 

4 . 38* 10~ 5 

13.90 

4.S3 x 10" s 

, 15.90 

^ -» J 

5. 35 x 10~ s 


T able' 2-b Entrainment' mas's flow per’ radian' (^')' at the 
outer ( BC ) boundary . 


* This was .confirmed by tests in which the mass flows 

were changed to three times the values tabulated. 
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At the upper boundary (where outflow occurs), and at the 
symmetry axis, the lateral velocity (v) is zero, and the 
inwardly-directed flux of u, h, and t is zero. 

Finally, the boundary conditions for R x and R^ at the 
symmetry axis are deduced from the knowledge that the net 
radiative heat transfer across the axis of symmetry is 
zero; elsewhere the presumption is that outgoing radiation 
escapes without reflection, while the incoming radiation 
is aT^ , where is the temperature of the fluid at the 
boundary. 

2,2 Auxiliary Inputs 

2.2.1 Thermodynamic' Functions 

The enthalpy (h) is calculated as follows : 


h = Zh-m. + H m fu 

( 2 . 2 - 1 ) 

where H is the lower value of the heat of formation of 
methane; the value used here is H=4.10 7 J/kg. The 
summation sign denotes summation over all species present. 
The species enthalpies hj are defined as follows : 


L. = ( C to 4- qj +c, 1 T l + c i3 T 3 )T 


( 2 . 2 - 2 ) 


The coefficients C- are tabulated in Table 2-C. 

- J 



18 



bhb 


mm 



(J/kg) 

(J/kg°K) 

(J/kg °K=) 

(J/kg°K 9 ) 

N 2 

10.435 

-5.3981 

1.2784 

-3.3472 

°2 

. 8.7061 

15.800 

-0.21869 

-0.13947 

o 

o 

to 

6.4748 

70.181 

-3.2833 

5.8576 

HgO 

18.343 

0.019232 

2.7700 

-7.2523 

(Vapour) 

ch 4 

15.0995 

212.57 

-1.8268 

-8.3016 


Table 2-C : The values of the Cj in Eqn. (2.2-2) 


The species concentrations in Eqn. (2.2-1) are deduced as 
described in Section 2.24 below, 

2.2.2 The Calculation of Viscosity 

The local mixture viscosity (p) is determined from the 
following formula : 


/* = fV + m c » 4 h,, 4 +(V+^+'V)p » 4 <2,2 " 3) 

Here, it will be observed, the combustion products are 
lumped together with the nitrogen; the error in ju. introduced 
by this practice will be small when, as is the case here, 
the nitrogen is present in larger quantities than HgO and COg. 
In any case, the viscosities of HgO and COg are not markedly 
different from that of Ng. , 
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The viscosities n ox and are related to the local 

9 

temperature via Sutherland's Law, which is : 


|*,= AjT 3/ 7^B j . + T) (2.2-4) 


The coefficients A- and B- are tabulated below : 

■J J 


Species 

Aj x 10 6 

— - - ^ q "Vi 

(kg/msec K ) 

B i <° K > 

■ CH 4 

1.00490 

171.473 

°2 

1.74410 

144.207 

N 2 

' 

1.40510 

111.452 


Table 2-d : The coefficients in Sutherland's viscosity 
law - Eqn, (2.2-4) ■ 

2.2.3 The Calculation of Density 

The local gas density (p) is calculated from the perfect- 
gas law : 


p = f/{R T S( m VM,)} 
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where p is the mean fluid pressure, presumed for the present 
calculation to he that of the ambient surroundings; R is the 
universal gas constant; and Mj is the molecular weight of 
species j. The summation is over all species present. 

2.2.4 The Chemical-reaction Models 


The chemical reaction considered is:- 


C H 4 + 1 0 2 * 2H,0 + CO, 


The reaction is presumed to he irreversible, and to occur 
in a single step. 

Two reaction models have been employed. The first presumes 
that the rate of reaction is infinite: fuel and oxygen 
therefore cannot coexist. When this is so, there is no need 
to solve the transport equation for m. ; both m, and m 

X Li 1 U Ua 

may be deduced from the mixture fraction f as follows* : 

O « f < f £t ; 


m fa = ° » Tn-o*,ou(l- j-') (2.2-6) 

1st ' 

if f st 4 £ 4 1 

O , ™£u = £ I (2.2-7) 

i - ist 

Here, m Qx is the mass fraction of oxygen in the ambient air, 
and f st > the mixture fraction for a stoichiometric mixture, is 


* These formulae, and others given later in this 

section, presume that at the jet is unity; i.e. 
■> the jet is composed of pure methane. 
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f = l/( 1 + i/tn ^ (2.2-8) 

st / \ ' ox, co y 

where, i is the mass of oxygen needed just to burn unit 
mass of fuel to completion; for the reaction considered here, 
therefore, i = 4. 


The above formulae are depicted graphically in Tig. 2-b, 


which shows the variation of 
and M 


“fu’ 


m ox ( an£ * also m. 


N, 


-h 2 o 



Fig. _2-b 


Variation of species concentrations with f for 
an inf initely-f ast reaction. 


In the second reaction model, the -rate of reaction is 
calculated from an Arrhenius expression. .The transport 
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equation for ra^ given in Section 2.1 is solved, and nt QX 
is computed from : 


m ox 0 £ ") - \[l m £ J (2 .2-9) 


Por both reaction models, the following formulae are 
employed for calculating the remaining species concentrations : 


^0 = 2 ( 2 . 2 - 10 ) 


^CO z — (f ~ (2.2-11) 


m w = 1 - TT~L _ TTL _TTl _ TO. , 0 , 10 . 

Ng. fu ox co i \o (2.2-12) 


2.3 The Solution Method 


2.3.1 The Finite-Difference' Equations 


The equations are solved by a marching-integration finite- 
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difference method, based on the work of Patankar and 
Spalding (Ref. 9), in which the domain of solution in 
the x-r plane is imagined to be covered by a mesh of 
grid points as shown below* : 



1 

■ 

■ 

■ 

■ 

a 

■ 

a 



8 

■ 

■ 

■ 

m 

a 

a 

a 



■ 

8 

■ 

a 

a 

8 

a 




■ 

a 

B 

B 

8| 

B 

a 

■ 




1 

a 

■ 

a 

B 

a 

a 





■ 

■ 

a 

a 

a 

a 



■ 


■ 

■ 

a 

a 

a 

a 



x 


Fig. 2-c : The finite-difference grid. 

The value of rf> at grid point P and at a particular time 
t n is related to the value at P at the previous time step 
(4>p), and to the values at the' neighbouring grid points 
X+, X-, Y+ and Y- , by integrating the general <j) equation 
(Eqn. (2.1-3)) over the small control volume or cell 
surrounding P shown in Fig. 2-c. For performance of the 
integration, realistic assumptions are made regarding the 
variations of in the x, y and t directions between grid 
nodes; for most purposes the variations are taken to be 
stepwise. The convective and diffusive fluxes are, however, 
calculated using a special hybrid difference scheme, due to 
Spalding (Ref. 10), which is designed to improve accuracy 


Because the general solution method described 
here is valid for Cartesian (x,y) coordinates as 
well as the cylindrical coordinates used here, 
J the symbol y is often used throughout this 
section to refer to the radial direction. ' 
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and to eliminate instabilities which can occur when the 
convection terms are large compared with the diffusive 
ones. The 3(($)/8t term is evaluated on the presumption of 
linear variation of $ with t; for all other purposes, the 
new values of <f> are assumed to prevail over the whole of 
time range t° to t n . 

The outcome of the integration is a finite-difference 
equation of the form : 

Ap = A. x+ <j> x+ + A x _4? x _ + Asy + 4> y+ 

+ + Ap 4>° +• S U ( 2 . 3 - 1 ) 

in which equation the coefficients A x+ , A x _, A y+ and A y _ 
contain the effects of convection and diffusion, while 
accounts for the unsteady term. The integral source 
te-rui has been written in the linearised form S u + S p^p*' a nd 
the coefficient S has been combined into A , which is 
defined as : 


A P - A X* + A *- +Ay + t A v _ + A° p ~'Sp (2.3-2) 

The velocities are treated, in two respects, differently 
from other variables. Firstly, the velocity grid points 
are located midway between the nodes for other variables. 
This arrangement, which is illustrated in Fig. 2-d, has the 
convenient feature that the velocities are stored- midway 
between the pressures which drive them. The second special 
feature of the velocity equations, is that, in preparation 1 
for the calculation of pressure by the means described in 


^ Footnote : The use of this form, when S is negative, is often 
found to enhance the stability of the solution; for, when (S ) 
is large, equation (2.3-1 ) tends to : p 

*pi>p + (“Vi = A P+r 


r Po 
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the next section, the pressure-gradient part of the source 
term is separated from the rest. The finite-difference 
equation for u, for example, therefore has the form : 


Ap Lip = _)_ A;(_ U x _ -j- A y+. Uy + 4- Ay_ Uy_ . 

+ A p U P + Sli -j- Dp pp} (2.3-3) 


The subscript nomenclature adopted for the velocities is 
illustrated in Fig. (2.3-3). Thus, P refers to variables 
stored at nodes within the boomerang-shaped envelope in the 
figure. 

Finite-difference equations may be derived for R and R , 

■ x y 

in a similar manner to that for 4 by integration of Eqns. 
(2. 1-8) and (2.1-9). 

The resulting equations have the form : 

— A X4 . Rjc, X_|_ "T" A x _ R.x_,)(- S u. 

A p Rjp'= Ay + R 3 ,y + 4- Ay„Ry,y_ 4- Su 

2.3.2 ’ The calculation of pressure 

Pressure is calculated by a novel procedure due to Patankar 
and Spalding (Ref. 9), which at each time-step the momentum 
equations are first solved for an estimated pressure field, 
which is subsequently corrected so as to ensure that the 
corrected velocity field satisfies mass continuity for each 
cell. This is achieved as described below. 


( J-3 - 1) 

( 1 . 3 - 5 ) 


The finite-difference form of the mass-.conservation equation 
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may be written as : 


LTe — fjO SxS^. + a K+ o u 

( t n - t° ) L 

+ . t V(V v * ~ ‘VfV'^P 


«p 

o (2.3-4) 


where Sx and <$y are all dimensions, and a , a ^ and a 

- x y+ y— 

refer to cell face areas. All these quantities are shown 
in Fig. 2-c. 


The velocity corrections (u' and v') can be related by the 
following approximate linear relations to the pressure 
corrections (p' ) : 


u' P » e ;c p;. - Pp) 

iTp = E p ( py_ - p; ) 


In steady-state calculations, the coefficients E^ and 

are simply set equal to D^/A and Til / A . For transient 

r r r r* r 

calculations, however, a stability-enhancing form* is used, 
in which , for example, is set to D^/Ap. 


The steady-state version, in contrast, presumes that 




» u\ 


= u’ 


y+ 


U ’y- = 0- 


Substitution of these expressions (and similar formulae for 
u’ x+ and u’y + into Eqn . (2.3-4), yields an equation for the 
pressure correction p’ , of the following form : 


*Footnote 


U 


/ 


: Which results from presuming that in a transient 
flow in which changes are taking place 
systematically, it is reasonable to suppose that 






w 'x- & U W 
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Ap - A x+ + A x - A A y+ f^y + -f- A v _fjy_ _ hip 

(2.3-7) 


where m is the mass source at P associated with the 
velocity field (u* and v*) calculated from the estimated 
pressures; i.e. : 

“ itl — L ? 2 + djt Px-Up 

(t n - t B ) ' 1 

.+ a V - <**-(%- (2.3-8) 

2.3.3 Solution of the Equations 

The overall— solution scheme is as follows; at each time 
step, the following operations are performed 

(i) The first operation is to calculate the fluid 
properties y and p at each grid node, from the 

temperatures at the previous time step. 

(ii) The momentum equations are then solved for the 
estimated pressure (p*) field, to yield u* 

and v*. The p*'s are usually taken as equal to 

o 
P • 

(iii) Next, the p 1 equations are solved, and the 
appropriate corrections are applied to the 
pressures and velocities, according to Eqns. 

(2.3-5) and (2.3-6). 

(iv) Finally, the equations are solved for the remaining 

<p’s. 
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In some cases (detailed in Section 3), the practice was 
adopted of iterating on u, v and p by repeating steps 
(ii) and (iii) several times, and using, for all iterations 
but the first, the calculated pressures at the previous 
iteration as the estimated pressure field. 

Finally', it should be remarked that the sets of finite- 
difference equations of. the form of Eqn . (2.3-1) are 
solved by the fast alternating-direction TDMA method 
described in Ref. 9. 

2.3.4 The Provision of Steady-state Solutions.' 

The initial conditions of the zero-gravity transient 
calculations, were solutions obtained for the steady-state 
situation with gravity for the same conditions. For these 
steady-state calculations, it is of course desirable to 
employ realistic initial guesses to the fluid-property 
fields; in particular, high initial temperatures must be 
provided in the region of the flame', in order to ensure that 
ignition occurs. For the computations ■ reported here, the 
initial conditions for the steady-state solutions were 
calculated by presuming parabolic profiles of u and uniform 
profiles of all other quantities within the flame (Fig. 

2-f ) . The overall balance equations for mass, momentum, 
mixture fraction and enthalpy for the conical-shaped flame 
are then : 


? 5 R 2 = 
1 2 


x R X + 

■ Y . 

(2.3-1) 

f u 1 R 1 = 

4 




1 2 

3 

l i y 

o 

(2.3-2) 

fu- Fff = 
l' , 


\ l / 

1 

1 

O 

(2.3-3) 
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p. u R b - ( ? ^ -v- Y \o (2.3-4) 

2 V 2 'o 

Here, the overbar denotes mean values within the flame, 
and subscript o refers to conditions at the nozzle; ip 
is the entrainment mass flow given in Table 2-b. 

These equations are solved, assuming infinite reaction 
rate, to yield H, u, h and f , which provide initial values 
for u, h , and T. The pressure is presumed uniform, 
initially, at atmospheric pressure, and v is set so as 
to satisfy local continuity. 

Furthermore, when the mj u equations were solved, a steady- 
state solution for infinite reaction rate was used so as 
to provide initial conditions for the finite-reaction- 
rate calculations. 

2.4 The Computer Program 

The calculations have been performed with a version of the 
CHAM computer code, EASI (the name stands for elliptic 
axisymmetric integrator ) , which is a general computer program, 
written in Fortran IV, for calculating flow and heat and 
mass transfer in two dimensional flow with or without re- 
circulation. Both unsteady and steady flows may be handled, 
and the solution may be performed in cartesian or cylindrical- 
polar coordinates . 
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3. DETAILS OF THE COMPUTATIONS 


3.1. Introduction 

Task II of the contract required performance 'Of computations 
for 25 different cases to assess the capability and accur- 
acy of the solution procedure. Some of these cases were 
devised to illuminate the effects of certain parameters of 
the numerical procedure, such as the number of grid points, 
time step and under-relaxation factors. Other cases inves- 
tigated the response of the solution to variations in the 
experimental conditions or the reaction-rate constants. 

'In the course of carrying out the computations, and due 
to. the. lack of experimental data concerning the reaction- 
rate constants, it was decided to perform computations for 
some additional cases. In fact a total of 21 extra cases 
were • studied, bringing the total number of cases to 46. The 
specifications of all these cases are presented below. 

3.2. _ Specifications of the Cases Studied 

Case 0 

The experimental conditions of this case are those specif- 
ied in Task I of the contract and thus the label 'standard 
ease' will be used later in the report to refer to this 
case • 

Case i 

The number of grid nodes was decreased by 51% from that of 
Case 0. 

Case 2 

The number of grid nodes was increased by 40% from thatof Case 0 
Case 3 

The standard time-step (.01 sec.) was decreased by a factor 
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of 3 for the first 0.4 sec. of the transient solution. 

Case 4 

The standard time-step was decreased by a factor of 10 for 

the first 0,4 sec. of the transient solution. 

< 

Case 5 

The standard time-step was decreased by a factor of 50 for 
the time interval starting at ,01 sec. and ending at 0.308 
see. 

Case 6 

The standard time-step was increased by a factor of 10 
throughout the interval of zero to 1 sec. 

Case 7 

The position of the cylindrical boundary BC (Fig.1) was 
moved outwn_rd by 50% of the standard value. The converged 
steady-state ( 1-g ) solution for the standard case was 
read off disk and employed as the starting conditions for 
this case, except that the velocities along the cylindrical 
boundary were appropriately modified to ensure the same 
fluid entrainment for both this case and the standard case. 

Case 8 

The position of the cylindrical boundary, BC was moved 
inward by 50% of the standard case. 

Case 9 

The number of iterations, at each time step, on the hydro- 
dynamic variables u,v,p were doubled. Previously, only 
one iteration was performed during the transient solution. 

Case 10 

As for 9 except that the hydrodynamic variables were iter- 
ated upon four times. 
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Case 11 

This case examines the effect on the transient solution 
of under-relaxing the dependent variables. The following 
formula for under relaxation was used to calculate the 
value of a dependent variable 0^ : 

- 0 -*)<% ■ + «**> 

where 0 and 0° are the values of 0 at the current and 
P P 

previous time-steps respectively/ and cC, is the under-rel- 
-axation factor. The value of cL was taken as 0.5 for all 
variables. 

Case 12. 

As for case 11 except that the value of cC. was decreased to .25 
Case 1 3 

The order of solution of the finite-difference equations 
was altered as shown in table (3-b) . 

Case 14 

A further change of the solution order was made as shown 
in table (3-b) . 

Case 15. 

The fuel flow-rate was reduced to 1,08 cm /sec. The long- 
itudinal length, X^, of the integration domain was reduced 
to '4 cm to accommodate appropriately the shorter flame. This 
adjustment of was based. on the analytical solution of a 
constant-property laminar diffusion flame that predicts that 
the flame length is proportional to the flow rate. The hydro- 
dynamic variables U, V, P were iterated upon twice to smooth 
out the solution. 

Case 16 

As for case 15 except the fuel flow rate was increased to 
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3 

5.15 cm / sec and the value of was correspondingly 
adjusted. 

Case .17 

Both the mass flow-rate and the nozzle diameter were 
increased as shown in table (3^-a). The longitudinal length, 
X^, of the integration domain was estimated as discussed 
earlier^ The radial width, Y^, of the domain was also estim- 
ated from the constant-property analytical solution, which 
predicts that the maximum flame width is proportional to the 
nozzle diameter. 

The contract required that for this case, the converged 
-.solution of the steady (1-g) condition be obtained and used 
as the initial condition for the following unsteady case. 

Case 18 

The experimental conditions for this case were the same as of 
Ease 17 except that here the gravity was zero. 

Case 19 

The fuel flow-rate was smaller than that of Case 18 but 
the nozzle diameter was increased as given in table (3-a), 

Case 20 

The fuel flow-rate and the nozzle diameter were both smaller 
than those of case 19. 

Cases SI to 25 

In the previous 20 cases, the infinite-rate reaction model 
(section 2) was employed to obtain the local mass fractions 
of fuel, oxidant and products. This model, because it assumes 
that reaction will always take place whenever fuel and oxidant 
co-ex'ist locally, cannot predict the flame extinction which 
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occurred experimentally in Case 20. Thus, in order to 
predict this extinction, we must employ the finite-rate 
reaction model (Section 2. 24), provided that the approp- 
riate values of the constants A and p are known. Since 
there was no suitable information to assist in ascribing 
the values of A and p these quantities could be treated as 
parameters that could be freely varied until the desired 
solution was obtained. The converged (1-g) infinite-rate 
solution was taken as a starting point; this corresponds to 
A=0. Then successively increasing values of A were employed 
(Cases 22A to 221). in a search for the critical value of A which 
just, kept the (l-g)flame alight, The idea was that the transient 
solution would predict the flame extinction if the critical 
value of A was employed. Cases 24A to 241 were also performed 
in the search for the critical value of A within a wider 
range. Although this critical value of A was not obtained in 
the few runs which were made, valuable experience was gained 
and will be reported upon..; it establishes necessary found- 
ation for future use of the solution procedure. 

The values of A used in Cases 21 to 25 are listed in table 
( 3-C) . 

i 

The purpose of Case 21 was to demonstrate that the finite- 
rate model when employed with a sufficiently small value of 
A (i.e. very fast reaction) will produce results which are 
consistent with those of the infinite-rate model. 

Cases 23 and 25 are both >for the transient study of two 

4 4 

different values of A, namely 10 and 1.5x10 . 

It can be seen from table (3-C) that the temperature was 
under-relaxed in Cases 21, 24A to 241 and 25. This practice 
was necessary to dampen the oscillations of temperature as 

This information was received through a telephone call .from 
Dr. Cochran . 
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will be discussed later. 

Case 26 

In this case, the Prandtl number, was increased to 0.7 
to examine its effect on the solution by comparing the 
results with those of the standard case where was 0.6, 

After presenting the specifications of the studied cases, 
we can proceed now to discuss the results. 



, Case 

Experimental cond's. 

IBi 


Time-Step Information 

Floy; Rate 
(cm 3 /sec) 


mmm 

IBBI 

mm 

MBH 


0 

2,8 

,102 

8. 

4,5 

.01 

0. 

1. 

1 

tt 

II 

TT 

ft 

11 

II 

tt 

2 

tr 

ir 

U 

tt 

11 

tr 

IT 

3 

ft 

tt 

II 

ft 

.003 ) 

o. 

0.4) 






.01 ) 

.4 

1. ) 

4 

tl 

Tt 

II 

ft 

.001 ) 

0.0 

0.4) 






.01 ) 

0.4 

1. ) 

5 

Tl 

tt 

tt 

tt 

.01 ) 

0.0 

.01 ) 






.0002 ) 

o.oi 

. 308) 

6 

II 

tt 

tt 

tt 

0.1 

0.0 

1. 

7 

it 

11 

11 

6.75 

.01 

0.0 

1. 

8* 

tt 

11 

Tt 

2.25 

11 

11 

TT 

9 

it 

• ft 

II 

4.5 

tr 

ri 

11 

10 

tt 

1 ft 

II 

11 

tt 

tt 

11 

11 

tt 

Tt 

II 

It 

" 

ii 

11 

12 

Tt 

tt 

Tt 

II 

it 

Tl 

It 

13 

II 

11 

If 

II 

it 

tl 

It 

~f4 

tt 

tt 

tt 


it 

T1 

If 

15 

1.08 

tt 

4. 

tt 

ii 

It 

tr 

16 

5.15 

tt 

15. 

tt 

ii 

II 

tt 

17&18 

5.50 

.165 

15. 

7.5 

Tl 

tt 

ii 

19 

3,20 • 

.884 

12. 

40. 

tl 

TT 

tt 

20 

2.42 . 

.372 

8. 

20. 

11 

tt 

tt 

21 

2.8 

.102 

8. 

4.5 

.01 

0. 

1.03 

22A-I 

2.42 

.372 

S. 

20. 

.01 

0,0 

1. 

23 

tt 

- tt 

11 

tt 

tt 

tt 

II 

24A-M 

tt 

tt 

tt > 

11 

. tt • 

tt 

t» i 

25 

tt 

tt | 

tt 

11 

tt 

Tt 1 

tt 

26' 

2.8 

.102 

8. 

4.5 

.01 

0. 

1.03 


Error in Computation 
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Table (3-b) 



■ 

■ 

M 


Itera- 

Under-relax- 

Solution Order of 

Case 



time- 

ation factor 

^Dependent variables. 


■j 





for dependerr 








variables 


0 

17 

X 

25 

1 


1.0 

u,v,p,h,f ,R X ,R 

1 

12 

X 

18 

TT 


tl 

tt 

2 

24 

X 

25 

tr 


Tt 

TT 

3 

17 

X 

25 

it 


It 

Tt 

4 


tt 


TT 


Tt 

11 

5 


tl 


tt 


TT 

n 

6 


M 


Tt 


TT 

ii 

7 


rt 


11 


tt 

ii 

8 


tt 


T 1 


It 

IT 

9 


tt 


2 


11 

M 

10 


tr 


4 


11 

IT 

11 


it 


1 


0.5 

It 

12 


tt 


II 


0.25 

" 

13 


TT 


II 


1.0 

v,u,p,f ,h,R ,R x 

14 


TT 


" 


it 

R j h , R , f , v , u , p 







x y 

15 


Tt 


2 


M 

u,v,p,h, f,R x ,R y 

16 


It 


ft 


It 

tt 

17 


11 


TT 


Tt 

tt 

18 


11 


tt 


tl 

Tl 

19 


II 


tt 


It 

Tl 

20 


II 


II 


TT 

II 

21 


IT 


1 


TT 

II 

22A-X 






TT l 

It 








23 


ft 


1 


TT . 

it 



ft 




ft 

It 









25 


TT 


1 


TT 

11 

26 


ft 


tt 


TT 

11 


*In all cases, the grid lines were non-uniformly spaced. A geometrical 
progression -was employed in both the x-and y-directions. 
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Table (3-6) 


, CaseJ 

Arrhenius 
Constant A 
(J/Kg raol/°K) 

Number of 
Iterations 
at each storage 

Under-relax- 
ation, factoi 
for temperat 
ure 

" 

-■Transient 
Studied? ■ 

21 

10. 

80 

= 3 

Yes 

22 A 

10,000 

10 

= 1 

No 

B 

11,000 

It 

tt 


C 

12,000 • 

>. It 



D 

13,000 

tt 



E 

14,000 

t? 


tt 

F 

15,000 

M 


1! 

G 

16,000 

Tt 



H 

17,000 

tl 


It 

I 

18,000 

11 

II 

It 

23 

10,000 

50 

m 

Yes 

24 A 

10. 

100 

= 3 

No 

B 

5,000 

20 


t» 

C 

8,000 

tt 

TT 

1! 

D 

10,000 ‘ 

Tt 

TT 

, TT 

E 

11,000 

tt 

tt 

tt 

F 

12,000 

tt 

tt 

tt 

G 

13,000 

11 

II 

II 

H . 

14,000 

It 

II 

II 

I 

15,000 

Tt 

tt 

Tf 

J 

16,000 

tt 

tt 

tt 

K 

17,000 

IT 

> 


tt 

' L 

18,000 

It * 

»T 

tt 

M 

19,000 

tt 

tt 

T! 

25 

15,000 

100 

= 3 

Yes 


For cases 22 and 24, a stage refers to a new value of A, while 
it denotes the ^steady (1-g) condition for cases 21,23 and 25. 
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4. ANALYSIS AND DISCUSSION OF THE RESULTS 
4.1. Introduction 

In this chapter, the results of the computations are pres- 
ented and discussed by the aid of figures 4-a to 4-i. and by 
reference to the attached computer print-outs. 

All figures, except F-ig. 4-h, compare the variation with time 
in zero gravity of the flame length for the particular Cases 
studied. 

From considerations of dimensional analysis discussed in 
detail in Chapter 5, it was found that the transient beh- 
aviour of the flame can be best studied by expressing the 
flame length and time in a dimensionless form as and ? 
respectively. 

|V ^ 

The definitions of Xp L and t are: 



-y r — 7 


Here, is the flame length defined as the distance 

along the axis of the fuel nozzle where the mixture-fraction 
■f is equal to 

4.2, Results and Discussion 



The dependence of the solution on the number of grid points 
in the integration domain, is illustrated here by way of 

/V 7* , 

the variation of with t . This is done for the 
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three Cases 0, 1 and 2. The experimental conditions for 
these cases are identical. The finite-difference grids 
of the cases are 17 x 25, 12 x 18 and 24 x 25 respectively. 

Figure 4-a displays the transient behaviour of Xfi for the 
three grids. It is demonstrated that for a fixed time-step 
increasing the number of grid points in the domain of integ- 
ration (i.e. refining the grid) leads to larger oscillations 
of the predicted flame-length These oscillations may be 
linked with the numerical instability caused by the relativ- 
ely large magnitude of the convection and diffusion terms 
as compared to the transient in the finite-difference equation 
when the time-step is large. 


4, 2, 2, Effect of variation of the time step 

Four Cases ( 0, 3, 4 and 6 ) are chosen to demonstrate the 
dependence of the solution on the value of the time step. 
These cases have the same experimental conditions; and' the 
parameters that affect the numerical solution are identical— 
except for the time step. The time steps for the four cases 
are .01, .003, ,001 and .1 sec. The variations of X ft 
with t for these cases are compared in fig. 4-b. 

It is indicated in the figure that as we decrease the time 
step (Cases 3 and 4), the oscillations of the predicted XfL 
are dampened. The result of -Case 5 (time step = .0002 sec) 
is not plotted in this figure, since the variation of 
is almost .identical to that of Case 4. One may, however, 
compare the respecti%'e print-outs of Cases 4 and 5 after a 
time of 0,03 sec. It will then be clear that further ref- 
inements in the time step have insignificant effects on the 
solution. 
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The band of the experimental results is also shown in 
fig. 4-6. The predictions of cases 3 and 4 are seen to be in 
fair agreement. 

On the other hand, when the time step becomes very coarse(Case 6) 
large deviation from the meansurements is observed. 

4, 2.3. Effect of varying the lateral width of the integration domain 
The influence of increasing the lateral width, Y^, of the 
integration domain, is shown by Fig, 4-c, where the trans- 
lent behaviour of Xo is plotted for the two cases 7 and 
0. In Case 7, Y-^ is 50% greater than that of Case 0. It is 
illustrated in the figure that the transient variation of 
Xjrt, for Case 7 differs slightly from that of Case 0. 

We can compare the (1-g) flame-lengths from the respective 
print-outs of both cases. It is indicated that in Case 7, 
the flame length is 3.8 cm which is about 5,5% longer than 
that of Case 0 (3.6 cm). 

In the attempt of performing the computation for Case 8 

which examines the effect of decreasing Y^, a programming 
error occured in specifying the radial distance between the 
7th and 8th y-grid lines; the radial distance from the axis 
of symmetry of the 8th line was less than that of the 7th. 

This error produced an unrealistic result: the flame length 

after 1 sec in (Ovg) became about 63% longer than that of 
Case 0 (both cases have the same experimental conditions). 

The location of the error is marhed in the respective print- 
out . 

. It is expected, however, that when the grid is specified 
properly, a reduction in the value of Y^ will produce only 
slightly different results from that of Case 0. 
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4 . 2 . 4 . Ef feet ' of iterating at each time step 

rJ 

Here we compare the transient Xfi of three Cases (0, 9 
and 10) which differ only in the number of iterations per- 
formed at each time, step. Figure 4-d displays the results 
of the three cases. Two conclusions can be drawn.- 

First, it is shown in the figure, and confirmed by the res- 
pective print-outs, that increasing the number of iterations 
from 1 to 2 eliminates the oscillations in the transient 
flame-length. 

Secondly, it is clear that doubling the number of iterations 
from 2 to 4 produces identical results. • 

It is' interesting to note that the effect on. the transient 
Xft of iterating at each time step is almost the same 

as that of decreasing the time step (see figures 4-b and 

.... ■ * 

,4-d). It should be pointed out, however, that iterating 
twice at each time step is more economic than halving the 
time step since iteration applies only to the three dependent 
variables u, v and p whereas the finite-difference equations 
are solved for all variables at each time-step. 

4 . 2. 5, Effect of under-relaxation during the transient solution 

Three Cases (0, 11 and 12) are considered to study the effect 
of under-relaxing the dependent variables during the trans- 
lent mode of the solution. Again, the behaviour of Xpi is 
chosen to illustrate the sensitivity of the solution to 
the changes of the under-relaxation factor cC (see section 
3-2, Case 11). 

Figure 4 - e compares the transient for cases 0 and 11 
where the values of <x are 1 and .5 respectively. As 
compared to case O ( = 1), Case 11 shows that the time-rate 
of change of X// is much smaller, and gradually decreases as 
we march in time . 
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Of the three values oft* tested, the smallest (.25), case 12 
produces hardly any change in the flame length as time increases 
(see the print-out). Although it is expected that under- 
relaxation will slow down, to a certain extent, the progress 
of the redevelopment of the transient behaviour, the predictions 
show severe slowing down. The cause requires further 
investigation . 


4. 2. 6, Effect of changing the order of solving the equations 

Here we examine the results of changing the order of solv- 
ing the finite-difference equations of- the dependent var- 
iables. The behaviour of the transient flame-length for the 
Cases 0, 13 and 14 is displayed in Fig 4-f. The order of 
solving the equations in each case is given in table (3-b) 
and als.o in Fig. 4-f. 

The order of solution in Case 13 differs from that in Case 
0 in three respects: the equation of v is solved before that 

of u the equation of f is solved before that of -h, and, 

the equation of R is solved before that of R . Figure 4 exhib- 

, v X . 

its identical predictions for' the transient in the two cases. 

Now let us consider Case 14 and Case 0. In Case 14, the 
equations of R^, h, R^. and f are solved before those of v, 
u and pi whereas in Case 0, the equations of u, v and p 
were solved before the other four. The result of this diff- 
erence is displayed in Fig 4-f. It is shown that the tran?- - . 
ient flame-length in Case 14 is lagging, but otherwise iden- 
tical to, that of Case 0. This is' a consequence of the fact 
■ tha,t the flame length was determined (see Section 4.1) from 
the values of the mixture fraction, f,, which in Case 14 
was solved for prior to the hydrodynamic variables. The values 
of f, would be appropriate to the flow field (i.e, convection 
and diffusion) associated with the previous time-step. Thus, 
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for example , after the first time step' (.01 sec) of the 
transient mode, the flame length is that corresponding to 
the (1-g) steady state conditions. . 

We can then conclude that the solution algorithm is stable 
with respect to interchanges in the order of solving the 
finite-difference equations. 

1 

' 4. 2. 7. Effect of varying the flow conditions 

All the Cases (0 to 14) which have been examined up to this 
point, had the same flow conditions (i.e. the same nozzle 

O 

diameter, 0.102 cm, and s me methane flow-rate, 2.8 cm /sec), 

. No.w. the effect .of varying these conditions is studied by 
comparing the results, vs. t, of six Cases (15 to 20). 

The results of 3 Cases (15, 16 and 18) are plotted in Fig. 
(4-g) which also includes the result of Case 9. The experim- 
ental conditions of Case 9 were identical to that of Case 0, 
and in the computations of Case 9, the number of iterations 
(2) at each time step was the same as for the above three 
cases. 

Case 15 examines the effect of reducing the fuel velocity 

at the nozzle exit (U ) on the flame length. The flow rate 

of Case 15 is about 39% of that of Case 9, and both Cases 

_ ^ 

have the same nozzle diameter. Figure (4-g) shows that Xft, 
for Case 15 is about 2/3 that of Case 9. Furthermore, the 
value of the dimensional flame length, .0093 m, is in very 
good agreement with measurements (Ref. 7). At the end of 
the transient (0-g) mode, however, the predicted flame length 
is smaller than that for (1-g) condition; this result is in 
contradiction with experiment. More investigations are required 
to diagnose the sources of discrepancy in this relatively low 
Reynolds number Case. 

In Case 16, the effect of increasing U o on the flame length 
is studied. The flow rate of Case 15 is about .84% larger 
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than that of Case 9, and both Cases have the same nozzle 
diameter. Figure (4-g) demonstrates that, as expected, Xft 
at (1-g) condition is larger for Case 16 than for Case 9. 
Similar behaviour is also depicted throughout the transient 
mode. 

In Case 18, both the fuel flow rate and the nozzle diameter 
were respectively 96% and 62% larger than those of Case 9. 

The transient behaviour of Xft for the two Cases is displayed 
in Fig (4-g) . 

Two conclusions can be made from the results of Cases 9, 15, 
16,. 17 and 18: 


a) The predicted values ofx^for the range of flows studied are 
always less than $ (see Fig 4-g). This is in agreement with 
the results of dimensional analysis discussed in detail in 
Chapter 5. 

b) The predicted transient behaviour of Xm at (0-g) for all 

a/ 

cases indicates that up to a value of t of about 1.5? *fL 
decreases as t increases, then Xfu starts to increase until 
it reaches a maximum value. This Xf L value is, except for 
Case 15, larger than its value at (1-g) condition. 

The computations for Cases 19 and 20 met some difficulty. In 
these two cases, the nozzle diameter was respectively about. 

8.8 and 3.7 times, the diameter for the standard Case (O).This, 
as explained earlier in Chapter 3, required increasing apprecia- 
bly (10 and 5 times) the lateral width of the integration 
domain. At the outer cyclindrical boundary of this domain, 
the radial velocities were calculated from the values' of 
the entrained mass flow rates. These rates were based on the 
solution of GENMIX program (see Section 2.1). Although 
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for Cases 0 to 6 and 9 to 14 (nozzle diameter of .102 cm. 
and domain width 4,5 cm) the solution was not sensitive to 
the variation of these entrainment rates (see Section 2.1), 
it is evident from the print-outs of Cases 19 and 20 that 
the solution becomes sensitive to their values when is 
increased appreciably. This is demonstrated in the predicted 
distributions of the u-velocity at (1-g) condition of the 
two Cases. The respective print-outs show large negative 
u-velocities (recirculation) in the field, associated with 
a sudden acceleration along the axis and close to the nozzle 
exit. This, is attributed to insufficient entrained fluid at 
the cylindrical boundary. Detailed investigation of this 
problem is required. 

4 . 2 . 8. Examination of the chemical-reaction models 

The results obtained from the two chemical-reaction models 
presented in detail in Section 2.2.4, are discussed here by 
comparing Cases 0 and 21. The two Cases are identical in both 
t'he experimental conditions and the numerical aspects of the 
solution except for the reaction model employed. Case 0 
employs the infinite-rate model while in Case 21, the finite- 
rate model is used with small value (10,)of the Arrhenius 
constant to simulate a relatively fast reaction. 

The results of both cases are displayed in Fig. (4-h) where 
the variations. of with distance along the nozzle axis 
are plotted. The figure shows, as expected, that the rate 
of decay of m fu along the axis is faster in the infinite-rate 
model than in the finite-rate one. 

One can also compare the temperature distributions in the 
respective print-outs of the two cases. At a given point in 
the flow, the temperature predicted by the finite-rate model, 
as one would expect, is -less than or equal to that of the in- 
finite rate model. 
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The transient behaviour of the flame lengths in both cases 
are also virtually identical (see the print-outs). 

4, 2. 9. The flame extinction study 

The search for the reaction-rate constants that produce 
flame extinction must, of course, be accompanied by a suit- 
able criterion for extinction; one such, is that the maximum 
temperature in the flow field be less than a specified tem- 
perature (e.g. 1000°K). Although we were, in our search, un- 
able to obtain extinction, we can report a definite trend 
which indicates that a wider range of search for the reaction?- 
rate constants will produce the values suitable for the flame 
extinction. 

It should be pointed out here, that a problem was encoun- 
tered, during the extinction study, with the convergence 
of temperature. We can summarize the experience gained in 
treating this problem by stating that the slower the reaction- 

rate, the more under-relaxation is required to achieve 

convergence (see table 3-c for the values of oC ) . 

4 . 2 , 10 .Effect of ■ varying 

The two Cases 0 and 26 are compared to study the effect of 
% on the solution. In Case 26, has a value of 0.7 as 
compared to 0,6 for Case 0; the two cases are otherwise 
identical. 1 - 

The transient variation of Xpi, for both cases are plotted 
in Fig (4-i). Obviously, the values of the Pr number play an 
important role in the solution of laminar diffusion-flames. 
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GENERAL DISC.USSION AND CONCLUSIONS 


5 . 1 Comparison of Predicted Results with Expectations 
Based on Dimensional Analysis 

When a laminar steady jet is injected into an 
atmosphere at rest, the rate of flow of in jected-plus- 
entrained fluid at a section distant x from the 
nozzle is equal to Sryx, where y is the fluid viscosity 
(presumed uniform) and buoyancy effects are absent. 
[Note that this result,' suprising because it is quite 
independent of the rate of injection, is a consequence 
of boundary-layer theory; but this theory provides a 
good approximation in the Reynolds-number range in 
question here.] 

Nov/ 1kg of methane requires about 16kg of air for its 
complete combustion. We can therefore suppose that, 
approximately, the flame length will be given by 
equating Srry x.,-, to 16 times p u irD„ a /4; thus: 

co i, I JL O O O 



Actually the numerical value of this dimensionless 
flame length will differ from 1/2 because of:- 

(a) the action of buoyancy; 

,(b) viscosity increases within the flame; 

(c) incomplete mixing of the fluid at any section; 

(d) elliptic effects; 

(e) chemical-kinetic effects. 

Factors (a) and (b) tend to diminish the dimensionless 
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flame length, while factors (c), (d) and, probably, (e) 
tend to. lengthen it. The result is that (5.1-1) can be 
expected to give good order-of-magnitude estimates. 


The flame lengths predicted by the EGOMAC computer 
program are indeed of the order of magnitude indicated 
by equation (5.1-1); and they agree with expectations 
in showing that the dimensionless flame length is 
smaller, the greater the influence of buoyancy. The 
dimensionless quantity representing this influence is 
of course the Froude number which, for present 
purposes, can be conveniently defined as: 





Since we might as well define our typical height as 
the flame length given by the foregoing simple theory, 
we conclude: 



(5.1-3) 


The transient behaviour can be expected to depend 
. m 

upon a non-dimensional time t, for which the 

normalising quantity is the typical height, divided by 

the injection velocity; thus: 

, i 



(5.1-4) 


We therefore expect that the results for the transient 
flame length can be represented by way of the following 
relation between dimensionless groups: 


% = 0.1-5) 

with minor additional effects of Reynolds number, 
chemical kinetics, and changes in the temperature levels 
of the injected methane and the surrounding air. 
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The numerical predictions contained in the present 
report do confirm this expectation. 


5 . 2 The Influences of Grid and Other Features of the 
Numerical Solution Procedure 


When a phenomenon of this kind is computed numerically, 
the procedure itself introduces influences which can 
be expressed thus: 


Xfl ' f V; tjW, 





he ( 5 . 2-1 


This expression indicates that:- 

<v 

o the dimensionless time step, ot, influences the 
result; 

© the size of the domain of integration, represented 

' by r max/ D o and x max /x typical’ has soras efiects 5 
© the numbers of grid lines in each of the two 

co-ordinate directions, as well as departures from 
uniformity of spacing, also affect the predicted 
flame length; 

o the number of iterations per time step, and other 
details of the procedure such as the order in 
which the equations are solved at each time step, 
all have their influence on the computed results. 

Of course, with a satisfactory prediction procedure, fit!', 

N , N etc are all given such values as render their 
r x 

influence small; but what those values are can be 
established only by systematic exploration. 

The investigation reported in the present paper 
represents only an initial phase of the exploration. 
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However, it is evident that sufficient information has 
been accumulated to enable the next phase to be 
planned successfully. 


5 . 3 The Influences of the Properties Presumed 


If the grid were fine enough and the time step 
extremely small, we can be sure that the solution of 
the finite-difference equations would be identical with 
that of the differential equations; whether it would 
agree with experimental data is, however, an entirely 
different matter. 

In the EGOMAC computer program as currently constituted, 
the formulations of the thermodynamic, transport and 
■chemical-kinetic properties are simple ones, as is 
entirely appropriate to an exploratory investigation. 
Thus, the heat-conduction and material-diffusion 
processes are represented by constant Prandtl and 
Schmidt numbers : and the complexities of the real 
chemical-kinetic processes are all hidden within the 
two free constants of an Arrhenius reaction-rate 
expression . 

Because of this, and of the fact that there are many 
more parameters to vary even within this limited 
conceptual framework than could be investigated 
during the present explorations, it is not surprising 
that there are disagreements between the predictions 
made by the EGOMAC computer program and the available 
experimental data. Indeed, it is very gratifying that 
the agreement is so good, both qualitatively and 
quantitatively. 

What is now needed is a two-part investigation in which, 
first, the disposable constants of the present model 
are varied so as to optimise the agreement over the 
whole range of conditions, and, secondly, even better 
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agreement is sought through the introduction of 
refinements to the model. These refinements should 
include more realistic formulations for thermodynamic 
and transport properties and for reaction' kinetics . 


5 . 4 The Influences of the Boundary Conditions 

Problems of this "elliptic" kind require specification 
of boundary conditions at all boundaries of the domain; 
yet these are rarely known. 

In the EGOMAC computations, a fixed-rate-of-entrainment 
boundary condition has been imposed, the entrainment 
ratfe being chosen as that which a simplified analysis 
dictates. Although there is no reason to suppose 
that the nature of this boundary condition significantly 
affects the computation in the central physically- 
interesting region, a further investigation is needed 
before the point can be established with certainty. 


5 . 5 Computational Economy 

Because of the low limit on available funds for this 
project, no program development was done. The EGOMAC 
therefore does not contain all the devices for 
economising computer time and storage which are 
available in more recent CHAM codes, and in others 
under. development. 

It should therefore be mentioned that the EGOMAC code 
is rather expensive to run. For a given numerical 
accuracy, improvements could be incorporated which 
would reduce its running costs to below 50% of the 
present level. 
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5 . 6 Recommendations for Further Work 

This investigation has firmly established that 
‘computations of the kind performed by EGOMAC are 
potentially capable of correctly predicting the 
transient buoyant laminar diffusion- flame process that 
formed the centrepiece of the project; and of course 
it confirms belief that many other flame and flow 
processes can be predicted by extensions of the method. 
To this extent, the investigation can he counted as 
successful. 

However, the success requires to be exploited; and, of 
the various directions in which future work can go, 
the following appear to be especially fruitful 

(a) Improvements of computational economy 

Because it is now desirable to exercise the EGOMAC code 
extensively, the cost of doing so should be reduced by 
measures directed to improving computational efficiency. 
There are. many of these available; and they will become 
especially valuable when finer grids are employed in 
the interests of improved accuracy. 

( b ) Optimisation of physical inputs 

If a closer fit is required between the predictions and 
measurements for the transient methane-air flame, it 
will be necessary to optimise the various inputs. For 
example, the assumption of a single uniform Schmidt 
number for diffusion of all components is certainly not 
correct; it could be replaced by assuming individual 
'Schmidt numbers for each component; and of course the 
Stef an-Maxwell equations could be introduced as a more 
realistic alternative. 

If extinction phenomena are to be realistically 
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predicted over a range of experimental conditions, 
there is’ no doubt that an appreciably more complex 
chemical-kinetic model will be needed than the present 
one. Although it would have been possible, had more 
funds and time been available, to extend the search 
■for optimum constants for the Arrhenius equation, this 
would be hard to justify until the physical properties 
have been optimised. 

(c) Employment of simpler computational schemes 

It is possible, with hindsight, ,to recognize that the 
EGOMAC computer code is too complex and expensive for 
some of the tasks of interest to the sponsor. Thus, 
since /the parabolic equations represent a good 
approximation to the steady-state situation, and since 
such equations can be solved .(eg by way of GENMIX, Ref 
8) more expeditiously than the elliptic ones, it would 
make sense to optimise the transport -property inpbts by 
an extensive investigation of steady-state flames. 

(d) Application to different systems 

Now that EGOMAC has shown its potential for the 
methaneuair flame, it can evidently be employed with 
some confidence for other flame systems also. Thus: 

o The reactants can be changed. 

o The boundary conditions may be altered. 

q Gravity may vary with time in accordance with a 
different program from the all-or-nothing 
prescription of the present example. 

o Turbulence may be introduced. 

• ^ 
o NOX reactions may be studied. 

c Etc. 

Further, it should not be forgotten that EGOMAC is just 
one member of a family of computer codes, which includes 
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those capable of handling three-dimensional unsteady 
flames. 


5 , 7 Summary of Conclusions 

(a) The EG011AC computer code has been shown to be 
capable of meeting the sponsor's requirements. 

(b) Many more computer runs have been made than were 
required' by the client; and their results are 
reported and summarised in this document. 

(c) It has been shown that the experimental behaviour 
can be reproduced with fair accuracy, provided 
that the time step is sufficiently short. 
Iteration at each time step also improves 
accuracy. 

(d) Further improvement of agreement between 
predictions and experiment can also be achieved 
by modification of physical-property data within 
reasonable limits. 

(e) Exploitation of the success requires further work 
in respect of;- 

e Improvements to computational efficiency. 

© Optimisation of physical inputs. 

© Conduct of supporting investigations by way 
of simpler computer codes. 

Extension to different and more complex flame ' 
systems. 


© 
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NOMENCLATURE 


Symbol 


Meaning 


V, a y +> a y. 


A j’ B j 


A p’ A x+’ A x~’ A v- 

C j ,n 


D U , D v 


E 

E u , E v 


f 

B 

h 

H 


i 

I, J 


K, L 


radiation absorption coeff- 
icient 

cell-face areas 

constant in the Arrhenius 
reaction-rate expression 

coefficients in Sutherland 
law for viscosity of species 
j 

coefficients in the finite- 
difference equations 

coefficients in the express- 
ion for enthalpy of species 

j 

nozzle diameter 

pressure-term coefficients 
in the finite-difference eq- 
uations for u and v 

black body emissive power 

coefficients in the correction 
formulae for u and v 

mixture fraction 

gravitational acceleration 

enthalpy 

lower value of the heat of 
formation of methane 

' stoichiometric ratio 

•radiative heat fluxes in the 
positive and negative r dir- 
ection 

radiative heat fluxes in the 
positive and negative x dir- 
ection 


m 

M 


mass' fraction 
molecular weight 
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P 


M_ 

P 

r 

R 

R x* 

(K. 


R 


S 


u' 


S 

P 


t 

T 

u 

v 



oc 

r 

SX,SV 

p 

p 

^r 

4> 

TW 


mass source at grid node P 

static pressure 
radial distance 
flame radius 

radiation flux quantities 

universal gas constant 

scattering coefficient for 
radiation 

source term 

coefficients in the linear- 
ised source term 

time 

temperature 
axial velocity 
radial velocity 
axial distance 
flame length 
radial distance 

constant in the Arrhenius 
reaction rate expression 

under-relaxation factor 
diffusion coefficient 
cell dimensions 
laminar viscosity 
density 

Prandtl/Schmidt number 
Stefan - Boltzmann constant 
general dependent variable 

mass flow rate of entrained 
fluid along cylindrical boundary 
up to x . .. • 
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Subscript 


Pertaining to. 


co 2 

£u 

h 2 ° 

3 

N 2 

ox 

P 

St 

X+, X-, ,Y+„ Y- 


J. 

0 

DO 


carbon dioxide 
fuel (methane-) 
water vapor 
species j 
nitrogen 
oxygen 
grid node P 

stoichiometric conditions 
neighhourning grid nodes 
variable (3 . 

conditions at the nozzle 
ambient conditions 


Superscr j pts 


' Pertaining to 


n 

o 

«i V 


conditions at the present 
time step 

conditions at the previou; 
time step 

velocities u and v 


i 


rd 


corrections to estimated 
pressures and velocities 

estimated quantities 

mean values within the flame 

dimensionless quantity 
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10 20 30 40 50 % 


Pig 4~g: The effect of Nozzle Diameter and Volumetric 

Plow Rate 


* Arrows indicate x fl values at (1-g) condition. 



Mass Traction of Fuel on' tho axis 
at t - 0,03 sees 
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Axial distance (crcs) 


Fig 4-h 


A comparison be tween., a Fast Finite-Rate Reaction and a Infinite-Rate 


Reaction for the Flow Conditio* of Case O 
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APPENDIX 
THE PRINT-OUT FORMAT 


This appendix is intended as an aid to the reader of the 
computer output which accompanies this report. At each' time 
step or after each iteration, several quantities are print- 
ed, and these are now to be individually discussed. 

TIME - The period during which the flame has been subjected 
to zero-gravity conditions. In some of the cases, e.g. 
Case 0, iteration is at first employed to obtain a 
(1-g) converged steady-state solution and for this _ 
situation TIM = 0. In other cases a converged steady- 
state solution is stored ( as a permanent file on 
computer disc storage), and is used to provide initial 
conditions. 

LFLAME-The flame length along the axis of symmetry. The 
flame envelope, or surface, is defined to be the 
locus of points, in the flow field, for which the 
mixture fraction . obtains the stoichoimetric value. 
The^Tlame length is determined by interpolation 
between appropriate grid nodes at the axis of 
symmetry. 

WFLAME-The maximum width of the flame, This ' quantity is 

determined by finding the width of the flame envel- 
ope along all the grid lines parallel to the y-axis. 
For any given grid line the envelope width is obtained 
•by interpolation between the appropriate grid nodes. 
The maximum of these envelope widths is set to 
WFLAME . 

U AT LF-The axial velocity at the tip of the flame. This 
quantity is also determined by interpolation. 



The mixture fraction is printed out at 8 locations along 
the axis. These locations are taken as fractions and 
multiples of the flame length (e.g, at x = 0,5 x LF) 

The remaining quantity which is printed out at each time 
step is the temperature at the nodal points along the 
axis and along the grid line which corresponds to the 
maximum width of the flame envelope. 

Finally more comprehensive print outs are given at the 
following instants. 

t = 0.0, 0,03, 0.1, 0.3, 0.5, 0.7, 1.0 secs. 

The variables printed out are:- 

u, h, m fu , m Qx , T, R x , R y ,pp, each at approximately 

100 points in the flow field. 
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